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DOCUMENTATION FOR PROGRAM SIMICE 


Introduction 


Program SIMICE calculates the conditions under which a model should be 
tested to maintain similitude between it and a prototype. It attempts to 
maintain the following constant between the model and prototype: 


(1) collection efficiency, Em, 

(2) heat transfer flux, q/A, and 

(3) the flux of, liquid water approaching the body under test, 
grams/meter /sec. ; 


The program solves for the liquid water content (LWC), the free stream 
velocity (V_), and average droplet diameter (D) of the model from various 
input data. The equations modelling these phenomena are nonlinear in form 
and are solved by the well-known Broyden technique. 


The user has flexibility in the data input and output. The inputs from 
a file are checked for reasonableness. The user must input initial values 
of the independent variables which are at least reasonable approximations to 
the final solution or the Broyden procedure may fail. 


The program writes out the following: 


(1) intermediate results from the BROYDEN program to the file DELBRO.DAT; 

(2) any error messages or flags which may be encountered to the file 
SIMICY.ERR; and 

(3) the input data and calculated values to the terminal or lineprinter. 


This method of handling the output permits the user to obtain the principal 
output in a highly-presentable form while, at the same time, having available 
intermediate and other output for more detailed inspection. 


The program is highly structured, and cop.ously supplied with COMMENT 
statements and, therefore, easily understood. It consists of the main program 
SIMICE and 14 subprograms. Program SIMICE 


a) reads in the necessary starting data (Sie DATA INPUT section); 

b) sets up the system of equations to be sc ved; 

c) calls on subroutine DELBRO (the Broyden soutine) which solves the 
equations; and 

d) writes out the final results. 


The subprograms (described in greater detail below) perform various chores 
such as 


a) calculation of physical properties; 
b) calculation of values of such quantities as the heat transfer 
coefficient, hos collection efficiently, a etc. 


This structure makes it possible for the user to easily replace current 
correlations with improved versions if and when they become available. 
AIR FORCE OFFICE OF SCIENTIFIC RESEARCH (AFSC) 
NOTICE OF TRANSMITTAL TO DDC 
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Mathematical Background 


The basic procedure used to find the conditions for the model tests 


scaled-down from the prototype test results according to the principals of 
similitude is as follows: 


1) write three equations governing the relevant physical phenomena and 
their scaleY¥-down and 


2) solve these equations for the three unknowns. 

The three unknowns referred to are: i 
LWC = Liquid Water Content for the model. 
Vor = Free stream velocity for the model. 
De = Average droplet diameter for the model. 


Other variables are brought in as parameters which can be changed as 
necessary by the user of the program. 


The physical phenomena whose equalities are to be maintained between the 
model and the prototype are: 


(1) the collection efficiency, ED 
(2) the heat transfer flux, q/A; and 
(3) the flux of liquid water approaching the body under test. 


The equations which model these phenomena are nonlinear in form. They are \ 
solved by the Broyden technique (2) using a program developed by Dr. J.D. Perkins. 


To facilitate the numerical work the program referred to above uses the 


ratios 
X(1) = (LWC) | / (LWC) - 
X(2) = De 
X(3) = Vou! ot 


where the subscript f refers to the model and m to the prototype. 


The three equations are developed in greater detail below. 


Collection Efficiency = En The collection efficiency En is a function of K (9° 


En = £,(K)) 7 
where g 2 A/LK jeoession FOr Y 
ro) 8 De I 
= tis : \ § 
APM, =o SRe,) pric TAB ie 
and ysrqnnounced Ss nee 
gene pastafieathon——— 
n=2 Pwo a ee 
9 A. a eee ton / 
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We want (E,), = (E,) at nipity Codes 
m eee: rene 1 and/or 
a ie ancink 
(Ee will equal ED. if (Ko) (Kon *\oast S¥ 


4 In the equations below the subscript f(referring to the model)became P, 

So, K 

_ _OMg), | : | 

Og) mn 

or Xn SA, 

Ky O76)» 

: (Np op) 

or Vm Dm Lep £, Roop’ m = 


2 
Vp Dp Lem £2(NRoD) p 


Let Log/Lep * R and the above equation becomes 


_f£QRe)m _ 
“£(MRe)p 


or 


Heat Transfer Flux = q/A . After Messinger(6), the heat flux to or from the 
body under study (t, < 32°F) : 


q/A = h[(t_-t_ - FW +) + 2.90L, (Psitte )] | 
SER Oe ~ oe JC P : 

c pa oO 

i 

(C Gtet joins te) = ee) : 
oat. po me aE | 
We want | 


(q/A), = (q/A) 


or (q/a),, / G@/A), = 1 | 
1 


F, = 1.0 - (q/A) / (a/A) , = 0 


Water Flux = (LWC) (V,) 


We want 
(LWC) (V) = (LWC) (VD 


or 


The above equations appear as follows in the program: 


QAM = (q/A) 


RAM = (A/A.) 
F(1) = 1.0 - QaM/(q/A) | 
F(2) = R- X(3) *K (2) *#24RAM/ (AIA) 


£(3) = 1. - X(1)*X(3) 


Trial values of X(1), X(2) and X(3) are continued until these three functions are 
arbitrarily close to zero. 


SUBPROGRAMS 
Each subprogram is discussed briefly below. 
SUBROUTINES DELBRO, BOUNDR, and GJINV 


These three subroutines together solve systems of nonlinear equations 
iteratively. They are based on the Broyden technique (2) as modified by the 
More-Cosnard procedure(8). The program was written by Dr. J. D. Perkins, Im- 
perial College, London. : 


The user must supply the set of functions F(I), set up the initial values 
of the independent variables X(I), and input boundary values of _,» BOUND(I), 
within which the program attempts to find a solution. Other input parameters 
are explained by rather detailed comments in DELBRO. 


DELBRO may fail to find a solution for a variety of reasons: 
(1) Initial Jacobian approximation may be singular because of: 
(a) poor scaling of variables; 
(b) equations are not independent; 
(c) poor initial guesses. 
The program works best if the independent variables are kept between -1 


and +1 and the functions F(I) are kept near 1.0. Dr. Perkins (7) has suggested 
transformations of the type 


k 
= (X,-a)/b 


i 
k 
fi = f fc 


(2) Bounds may impede progress towards a solution. It may be that a solution 
does not exist within the bounds, or, a poor starting point may have 
been chosen. 


(3) The dimensions NW, NB are too small. See listing of program for re- 
quirements. 


(4) The program may not have converged to the precision specified (ERROR) 
within the maximum number of iterations specified (MAXIT). 


The subroutine BOUNDR checks the values of X, calculated against the boundary 
values supplied by the program, sets INFO = (Konp} equal to 5, writes out a warn- 
ing and returns control to DELBRO which, in turn, returns control to SIMICE. 


Subroutine GJINV inverts the preliminary Jacobian matrix. 


FUNCTIONS DENL, PSG, CW 


Function DENL calculates the density of liquid water (g/cc) from ~50°C to 
+50°C. It is based on data and equations given by KELL (4). Below -20°C the 
results are extrapolated. 


Function PS@ calculates the vapor pressure of ice (Kilopascals). It is 
based on the work of Arnold Wexler(13). 


Function CW calculates the heat capacity of liquid water in the range -50°C 
to +40°C. It ts based on data from Smithsonian Meteorological Tables (10). 


FUNCTIONS MUMIX, KMIX, NPRMIX 


Function MUMIX calculations the viscosity of air-water vapor mixtures. It 
is based on equations (9-3.9) and (9-5.4) of Reid, Prausnitz and Sherwood(9). 


Function KMIX calculates the thermal conductivity of air-water vapor mixtures 
from equations (10-4.1), (10-6.1), (10-6.2) and (10-6.3) of Reid, Prausnitz and 
Sherwood (9). 


Function NPRMIX calculates the Prandtl number of air-water vapor mixtures. 
It uses KMIX, MUMIX and heat capacity data on air from (5) and on water vapor 
from (3). 


FUNCTIONS LAMDAS, EM, EMC and EMS 


Function subprogram LAMDAS calculates a value of \4/d_ for a given value of 
the drop Reynolds number Npgp. The equations it uses are based on data tabulated 
in reference (11). 


Function EM calculates the collection efficiency, E,. It first calculates 
a value of K = 2/9pwVoD"/(ugL,), the Reynolds number and Ky = LAMDAS(Np,)*K. It 
then calls on EMC to calculate Em for a cylinder, EMS to calculate E, for a sphere, 
or EM@@12 to calculate Em for a NASA 0012 airfoil. 


Functions EMC, EMS and EM@@12 are based on equations obtained from curves 
given in references (1 , 11). 


This set of Programs should be replaced by a procedure which integrates the 
basic differential equations (12) to obtain the collection efficiency, En - 


FUNCTIONS HC, QA 


Function HC calculates the heat transfer coefficient he from equations given 
in (14). 


Function QA evaluates the heat flux as given by equations in Messinger (6). 


soe 


13. 


14. 


References 


Bowden, D.T., A.E. Gensemer and C.A. Skeen, Engineering Summary of Airframe 
Icing Technical Data, FAA ADS-4 (1964). 


Broyden, C.G., Math. of Computation, 19, 577-593 (1965). 
JANAF Thermochemical Tables (NBS), 2nd Edition (1970). 
Kell, G.S., J. Chem. Eng. Data, 12, 66-69 (1967). 


Kothandaraman, C.P. and S. Suframanyan, Heat and Mass Transfer Data Book, 
Wiley (1975). 


Messinger, B.L., J. of Aeronautical Sciences, 20 (#1), 29-42 (1953). 


Perkins, J.D., Private Communication, August 1980. 


Pfeifer, G.D. and G.P. Maier, Engineering Summary of Powerplant Icing Tech- 
nical Data, FAA-RD-77-76 (AD/A-045 087) (1977). 


Reid, R.C., J.M. Prausnitz and T.K. Sherwood, The Properties of Gases and 
Liquids, 3rd Edition, McGraw-Hill (1977). 


Smithsonian Meteorological Tables (1951). 


Sogin, H.H., A Design Manual For Thermal Anti-Icing Systems, WADC-TR 54-313 
(AD 90156) (1954). 


Tribus, Myron, G.B.W. Young, and L.M.K. Boelter, Trans. A.S.M.E., 70, 917- 
982 (1928). 


Wexler, Arnold, J. of Research of the National Bureau of Standards, 81A, 
5-20 (1977). 


Whitaker, S., Elementary Heat Transfer Analysis, Pergamon (1976). 


MESES Te ese ee 


DATA INPUT 


Data input is from a file named ICE*DAT. The input is unformatted: each 
number should be separated from others by commas or at least one space. Form 
and meaning of the input are shown in the following table: 


Record 


i 


me NM NR NY 


Variable Variable 
Name Type 
IB Integer 
IH Integer 

R Real 

LC Real 
ITPRINT Integer 
TS Real 
TO Real 
PMG Real 
vg Real 
OMEGA Real 
LWC Real 

DM Real 
LIM Integer 
TFS Real 
TFO Rea 
PPG Real 

x (Id. 

e (2). ---~>>Real 


Significance 


Body type = 1 for cylinder 
= 2 for sphere 
= 3 for 0012 NASA Airfoil 


Heating condition = 1 No heat input 
= 2 Heat input 


Scale-up Ratio = (LC) /LC) 


Prototype characteristic dimension, cm 


Variable controlling printout by 
DELBRO = -1 No printing 
= 0 Final results only 
= 1 X + F on every iteration 
= 2X + F + W on every iteration 


Prototype surface temperature, °C 
Prototype free stream temperature, °C 
Prototype free stream total pressure, Kp 
Prototype free stream velocity, meter/sec. 


Model(and prototype)specific humidity 
= (Kg water vapor)/(Kg moist air) 


Liquid Water Content, g/(cubic meter) 
Model droplet diameter, microns 
Program Control parameter: If 


LIM = -1 --- end of run 
# -1 --- expects another set of data 


Model surface temperature, °C 
Model free stream temperatue, °C 


Model free stream pressure, KP 


Initial values of ratio variables for 
BROYDEN 


Variable 


Record Name 


5 & 6 BOUND (1) 


BOUND (2) 


BOUND (6) 


Variable 
Type 


Real 


fe entice mM RMP ee 


Significance 


Boundary values for the X(I) variables: 
Lower bounds in BOUND(1), BOUND(2), BOUND(3); 
Upper bounds in BOUND(4), BOUND(5), BOUND(6) 


A sample set of input data is given on page 10 . 
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DATA OUTPUT 


Output from this program is illustrated on the following pages. The first 
set of output is the principal output which comes out to a terminal or line printer. 
The second set shows the intermediate results which can be written out to the file 
DELBRO.DAT. The third set shows various intermediate results which can be written 


out to the file SIMICY.ERR. 


RUM St Leer 


RUN NUMBER= 


Le = FPRoTOy 


PROGRAM SIMICE 


1 BROYDEN CONVERGED 


2B-JAN-31 


MONE! AND FROTOTYFPE BONY CY!INTER 


FROTOTYFECINFUT: CONDITIONS: 
YRE CHARACTERIGTIC UIMENSTON= 


TMS= FROTOTYPRE SURFACE SEMFERATURE= 


TMO= PROTOTYPE FREF STREAM TEMFERAGTURE= 
TYRE FREE STREAM JOTAL FRFSSURE= 


PMO= FRUSUT 


Yao FROSOT 


Lutes PRUTOT 
A = ROTO 


LCOF= MOTE 
TFS= MODEL 
Tre= MODEL 
FRO= MOUE}. 
VEOs MOTE 
LWCF= MOLFL 


DF HMQNEL WATER URAPLE i 


OME GAs 


KEEN NIUMRE fs: 


Cc = FROTOI 


TRE 


YRe CREE STREAM VEL OCTTY= 


TRE TWO WATER CONTE UT 


TRE UROPLET UIAMFPTER= 


AGEL. CONDT TANS & 
CHARACTERISTIC DIMENSTOrd= 
SURFACE TEAPERATURE= 
STREAM TEMFERATURE= 
STREAM TOTAL FRESSURE= 
FREPF STREAM VELOCTIVY= 

LiQuIb WATER CONTENT= 
CIAME FER= 


DAMWGE-OLCKG H2N VAPORI/ CRG 


~ BROYTEN CUNVERGEL 


AFTER 15 CaAlLs 
4.76CM 
0.00C 
-10.00C 


7. GOKF 
116. 79METER/SEC 

2. OSGRANS/CUBIC METER 
32.G3MICRONS 


2.00CH 
0.00C 
-16.00C 
4”. OOKF 
96. O1METER/SEC 
2,49GRAMS/“CURIC METER 
2O.GOMICRONE 


MOTST AIR: 


15 CALLS 


MODEL AND FROTOTYPER Bunyi cy. INGER 


PROTOVYPECTNF UTS CANT TIONS: 
CHARACTERIST IL biMeNSTON= 


TMS= PROTOTYPE SURF AE TRAPFRATUFE= 


TMO= FROTUT 
BMO= FROTO! 
VMNO= FROTOT 
LWP PRATOS 
Le = PRUTOT 


LOor= MOLE. 
TFS= MOTEL. 
YtO= MODEL 
FRO: MOREL 
VFO= MOTEL. 
LWCrF= MODEL 
DF=MODEL Wé 


JMEGA= 


FUN NUMBER = 


YRE FREE ST*FAM TEMFERATJRE= 
STREAM TUTAL FRESSURE= 
fF. FREG OTREAM VELOCITY + 

YEE Lraurtp WATER CONTENT= 


YEE FREE 


YRE [WROPLET DITAMETER= 


MONEL CONDITONS$ 


CHARACTERISTIC DIMENSTON= 
SURFACE TEMPERATURIZ= 
STREAM TEMPERATURE = 
STREAM TOTAL FPRESSURE= 
FREES STREAM VELOCTTY= 


LIQUID WATER CONTENT= 


TER UROPLET DLAMETER= 


3 BROYDEN CONVERGED 


4.76CM 

0,00C 
-160,00C 

7 0OKF 
120,00METER/SEC 

2-O5GRAMS/CUBIC METEK 
32. O03MICRONS 


2.00CM 
6.9000 

~10.00C 

47 .OOKF 

98.93METER/SEC 
2+49GRAMS/CUBIC METER 
19.97MICRONS 


GO 4C00E-O2(KE H20 VAFURI/ (KG MOIST ATR) 


AFTER 14 CALS 


MOnEL ANT PROTOTYRE BONY ICTL INGER 


PROATATYRE SC THEIUIT COATT T TANS $ 


LE =M CUE 
OMe ays 


RUN NUMBER? 4 


Lab = 
TR eis 
fh Os 
FRO: 
eo Oe 
WOE = 


Ve = MUTE L, 


reurud 
FROTOT 
PROTOT 


PROTO 


PROUT 


rReery 


MOTEL 
MOUEL 
MONTEL 
MOTEL 
MOTEL. 
MO CED. 


PROTOT 
PROTOT 


PROTO 
FRUToUs 
PROTA 


Pye? 


a 


MODEL 
MOTEL 
MUTED 
MOTIEL. 
MODEL 
MOTEL 


GMEGA= 


RUA 


NUMBE 


FRO POT 


EROTUT 


PHOTON 
PREG 
PRT CET 


PROTOS 


MOGEL 
MOLEL 
MOSEL 
ROME L 


ager tine 


WATER 


FREER STREAM 
YEE 
FREE E€TREAM VELOCITY= 
LIQUID WATER CONTENT= 
NWRUFEET UTAMETER: 


ee 


TEMPERATURE = 


vb te 
TRE 
YRE 
YE 
YP E 
YFE. 


CUNFT TONS « 


HOUFL. 


CHARACTERISTIC iIMENSION= 
SUFF ACE TEMPERATURE = 
STREAM TEMPERATURES 
STREAM TOTAL FRESSURE= 
FREE STREAM VELOCITY= 
"QUE WATER CONTENT= 


WATER URAFLET DIAMETER= 


OY, 20008 -O0°06 HE VAR ORD 


RROYREN COAVE RIE tt 
MOTE! ANU FRUIT OE TEE 
PROPOQPY RE CLARE G3 
YPE UCHARADTER LST EC 

TYRE SE GU TE MBE Re TURE = 
PFE FIRE 4ATNEAM 
YEE FREE STREAM FUTAL PF 
pee PREE SSF AM URLACT Ty = 

YEE TQ WATER CANTEEN Ts 

VRE CRORLE YS NTQMETFR= 


AMREL CONT TONS: 
UHARACTERISTIC DIRENS TON] 
SURFACE TEMPERA TUSE = 
STREAM TEAPL AR ATURE = 
Pe aM TOTAL RRP SSURE = 
OG PRR SELLE DT t= 
WATER CONTENT= 
ULAME TERS 


wit iy 
VWaAOELI & 


Oo. 40G0E-C2°0hG Hoth VaPGR1 


R= § HROYNEN CONVERGE N 
MODEL AND FROTOTYEE 
PROTOTYPE CTNFUT? 

YR TEMPERATURE = 

YrE 

VEE, 

YEE 

YRE 

Y FE 


SURE ACE. 
FREE 
FREE STREAM TOTAL 
FREE STREAM VELOCITY y= 
LEQuin WATER CONTEXT = 
TROPLE? UTAMETER= 


MOM CONTI TONS § 
CHARACTERISTIC DLIMENSTON= 
SURFACE TEMPERATURE = 
STRE GY TEMPERATURE = 
STREAM TO AL PRESSURE 


were AEA 2 8A bore ee hy 


StREAM TOTAL PRESSURE 


“FG 


PEMPERATURE 


PES URE = 


Pehle: 


STREAM TEMPERATURE = 
PRESSURE = 


wow 


~19.00C 

47. OOKF 

115. OOMETER/SEC 
2 OSGRANS/CURIC 
32.035MICRONS 
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METER 


2.00CM 
2600C 
-10,00C 
47. Q0KF 
94.37METER/SEL 
D2, SOGRAMS/CUBIC METER 
#0.Q02MICRONS 


MOLET ALR? 


ROMY SOY. INGER 
CONETITIGNS ¢ 
DUMENSITON= 


4.76CH 
4.900 

-19.00C 

47. 00KF 

119, 00METER SFL 
2. O5GRAMS/CURIC METE! 
50. COMICSONS 


2 OOCM 
0.060 
~10.00C 
47 ..00KF 
2ehITGRAMS. CUBIC METEL 
20 ,O7MICKUNS 
MOG? 


AIK} 


AFTER 13 CALL: 


BRODY s Crt INTER 
CEUNILT TONS § 
PROTOTYFE CHARACTERISTIC EIMENSTON= 


4.°76CM 

O.00C 
“10,000 

47 GOKP 
100.90METER/SEC 

= OSGRAMS “CURTC 
32. O3MICRONG 


MITER 


2.00CM 
C.Q0C 
-10,090C 

4/,00KF 


£IMOTI 


OA ose 


* 
Aa 
* 
* 
* 
* 
* 
* 
* 
* 


* 


2 + 


% * 
et 
% + 


44% 


+ 
% 
eH tht 


2 F 


* 
* 
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SEEK EKEKK CEKE 
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